Vector magnetic hysteresis of hard superconductors 
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Critical state problems which incorporate more than one 
component for the magnetization vector of hard superconduc- 
tors are investigated. The theory is based on the minimiza- 
tion of a cost functional C[H(x)] which weighs the changes of 
the magnetic field vector within the sample. We show that 
Bean's simplest prescription of choosing the correct sign for 
the critical current density J c in one dimensional problems is 
just a particular case of finding the components of the vector 
J c - Jc is determined by minimizing C under the constraint 
J 6 A(H,x), with A a bounded set. Upon the selection of 
different sets A we discuss existing crossed field measurements 
and predict new observable features. It is shown that a com- 
plex behavior in the magnetization curves may be controlled 
by a single external parameter, i.e.: the maximum value of 
the applied magnetic field H m . 



PACS number(s): 41.20.Gz,74.60.Jg, 74.60. Ge, 02.30.Xx 



I. INTRODUCTION 

The irreversible magnetization of type-II supercon- 
ductors was successfully treated by Bean's critical state 
model (CSM) in the early sixties.lil The basic idea is that 
the material develops a maximum (critical) current den- 
sity in those regions which have been affected by an elec- 
tric field, opposing to the magnetic flux changes. With 
remarkable simplicity, the CSM could explain the depen- 
dence of the magnetization on the macroscopic dimen- 
sions of the sample, as well as the observed hysteresis. 
Both features are a manifestation of the non-equilibrium 
thermodynamic processes which take place in the experi- 
ments and, thus, cannot be explained by the Abrikosov's 
flux line lattice (FLL) theory.El Furthermore, hard type-II 
materials develop such a pronounced hysteresis that the 
reversible contribution from the equilibrium FLL may be 
usually neglected. Nowadays, the hit of Bean's intuition 
is well understood in terms of the FLL dynamics in the 
presence of pinning centers. a In the process of punch- 
ing vortices in or out the superconductor one produces 
metastable equilibrium states for which the gradient in 
the density of vortices is maximum. This results in the 
development of the macroscopic critical current density 
J c and corresponds to the balance between a repulsive 



vortex-vortex interaction and attractive forces towards 
the pinning centers. 

Bean's ansatz has permitted the development of the 
theory with a maximum simplicity while incorporating 
the main experimental facts. However, it was early rec- 
ognized insufficient for problems in which lattices of non- 
parallel flux tubes must be considered. For such cases, 
the current density within the sample is no more a vector 
with fixed direction, and additional prescriptions are re- 
quired. On the other hand, a wealth of phenomena have 
been reported along the last decades associated to crossed 
field measurements. Refs. provide some outstand- 
ing experiments, and the interested reader is directed to 
the papers cited therein for more complete landscape of 
the topic. Roughly speaking, when the experimental sce- 
nario is such that compression and rotation of vortices 
are induced, both the scalar repulsion between rigid par- 
allel vortices and the cutting energy barrier for adjacent 
twisted flux lines come into play. In particular, this has 
consequences on the transport properties of type-II su- 
perconductors in a magnetic field, as the pinning .of the 
FLL is strongly influenced by its internal stiffnessfl 

Several efforts have been made for the interpretation 
of rotating field experiments. Bean himself developed an 
approach for such problems .13 However, to our knowledge, 
the most extensive generaUiritical state theory is due to 
Clem and Perez-GonzalezliTO who have developed the 
so-called double critical state model (DCSM). These au- 
thors have provided a physical basis for the limitations on 
the current flow both parallel and perpendicular to the lo- 
cal magnetic induction. Thus, the so-called parallel criti- 
cal current density J c \\ and perpendicular critical current 
density J c ± were introduced. J c » relates to the flux cut- 
ting threshold, whereas J c j_ stands for the conventional 
depinning current density. Recently, a more sophisti- 
cated approach was introduced, the so-called two-velocity 
hydrodynamic modelu Essentially, this theory incorpo- 
rates the flux pinning and cutting phenomena within the 
framework of FLL's which consist of two vortex subsys- 
tems. Under certain circumstances it is fully equivalent 
to the DCSM, but, different, non-expected scenarios are 
also predicted. ■— ■■— ■ 

In earlier worktUo we have envisioned the CSM realm 
as a phenomenological approach which can be formulated 
by means of a variational problem with constraints. This 
permits to deal with a wide class of restrictions on the. 
current density. Thus, the optimal control (OC) theoryEj 
allows the minimization of a cost functional C[H(x)] 
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which weighs the changes of the magnetic field vector 
within the sample under the constraint J € A(H,x), 
with A a bounded set. This general formulation gives 
freedom to fix the particular region A of feasible current 
densities. For instance, the isotropic hypothesis \J\ < J c 
corresponds to choosing A as a disk, and the DCSM con- 
ditions |Jj_| < J c ±, |Jii| < J c ii, to an oriented rectangle. 

In this work we present a detailed discussion of the 
variational principle (Secjlj]), as well as the general as- 
pects about the Hami lton ian equations which arise from 
our formulation (Sec. Scc.|y| gathers a series of sim- 
ulations of the field penetration profiles for two different 
selections of the control space A. In particular, it is ex- 
plicitly shown that the DCSM becomes a particular case 
of our treatment. Then, Sec. [v| is devoted to observe 
some properties of the magnetization loops which can be 
derived from the simulated profiles. In Sec. VI we give a 



general overview of the model and discuss the compatibil- 
ity of different choices of A with available experimental 
data. 



II. VARIATIONAL PRINCIPLE 

In this section we will show that the coarse-grained 
electrodynamics of type-II superconductors in the crit- 
ical state may be formulated as a minimization prob- 
lem. From a more conventional side, on using the DCSM 
approach, one can obtain the successive field penetra- 
tion profiles in a magnetization process by means of 
Maxwell differential equations, and the constitutive E(j) 
law for the material. These equations, together with the 
boundary conditions allow to solve for the dependencies 
H(x,t), J(x,t), E(x,t). As it is usual in hard supercon- 
ductivity we are assuming B = p^H. 

Now we recall that many physical theories, which 
are typically posed by a differential equation statement 
have been interpreted as the minimization of a given 
functional. For instance, the Lagrangian formulation 
of classical mechanics relies on the equivalence of New- 
ton laws and the stationarity condition for the action 
S = j Ldt = j(T — V)dt, when the time evolution of 
the system is determined. Below, we derive the func- 
tional C whose minimization is equivalent to the stan- 
dard approach for the critical state in superconductors. 
In order to gain physical insight, we will infer the CSM 
equations after considering some aspects of the more fa- 
miliar eddy-current problem in normal metals (E = pj) . 
Assuming a discretization scheme in which H n stands for 
the magnetic field intensity at the time layer n5t and 
Ampere's law (V x H — J), the successive field profiles 
in a magnetic diffusion process may be obtained by the 
finite-difference expression of Faraday's law 



(-ffn+i — H n ) 
Tt 



Vx£ = -pVx(Vx H a+1 ) , (1) 



which defines a differential equation for iT n +i in the 
spatial degrees of freedom. Then H n +i can be solved 
in terms of the previous field distribution H n and the 
boundary conditions at the time layer (n + l)St. The 
minimization statement can be found by inversion of the 
Euler-Lagrange equations with respect to some Lagrange 
density that belongs to Eq. (|l|) . The first term is straight- 
forwardly inverted by integration with respect to the field 
H n+1 . This gives n (Hl +1 /2 - H n ■ H n+1 )/St, i.e.: 



gCgg+i/g - H n ■ H n+1 ) 

dHn+l 



The second term is a bit more involved as it cannot 
be transformed to an exact derivative with respect to 
the fields. However, if one introduces the notation 
(V x H n+1 )i = (,,).(),- II,, . ..k (with e ljk the Levi-Civita 
tensor) it can be identified from 



da 



-2[V x (V x if n+1 )]i 
Thus, Eq.(Q) can be rewritten as: 



d(d Xj H n+lti ) 



(2) 



where 



•7"n+l — ^o(-H„ +1 /2 — Hn ' Hn+l) + 

P 8t{V x H n+1 ) ■ (V x i? n+1 )/2 . 

From variational calculus it is known that Eq.(^) cor- 
responds to the necessary first order conditions for the 
functional Cm = Jq^u+i to have a minimum with re- 
spect to the field H n+ i (Euler-Lagrange equations) 

Finally, it is apparent that for minimization purposes 
Cm may be completed by means of a constant term to 
the form 

Cm = I | -Hn+i ~ -H n | 2 + 

i j P 8t(V x H a+1 ) ■ (V x H a+1 ) , (3) 

as H n (x) is considered to be given. 

In conclusion, from the mathematical point of view, 
minimizing Cm for each advancing time layer is just an 
equivalent statement of Faraday's laws in time discretized 
form. In addition, Eq. (|^) allows a physical interpretation 
which can be used for comparison and generalization to 
other systems. Notice that the minimization of Cm bal- 
ances the screening term \H n+ i — H n \ 2 and the isothermal 
entropy production term (<S = E-J/T). This can be iden- 
tified as a quite general property of dynamical systems 
subjected to dissipative forces. The quasistationary time 
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evolution of the system holds a compensation between 
some inertia term and the irreversible loss of energy. The 
minimum principle for the global entropy production rate 
lies behind the previous statement. It was introduced by 
I. PrigogineEa in the context of linear irreversible ther- 
modynamics and applies for non-equilibrium stationary 
states. As a simple exercise, the reader can check that 
minimizing m\v n+ i — v n \ 2 + 7U n +i ■ v n+ i6t with respect 
to u n +i one reproduces the discretized statement of New- 
ton's law for the motion of a particle against a viscous 
damping force. In the case of eddy currents |i? n +i — H n \ 2 
stands for the magnetic field inertia, including both the 
magnetostatic energy term H 2 +1 and the term H n - H n+ i, 
which describes the work against electromotive forces. 

Eventually, we notice that hard type-II superconduc- 
tors may be treated by a modification of the functional 
Cm for normal metals. We want to emphasize that the 
concept of hard material is associated to a limiting case 
for the E(J) characteristic, i.e: the electric field is zero for 
current densities below a critical value (J c ) and abruptly 
raises to arbitrarily large values if J c is overrun. In fact, 
this vertical graph limit is attained for the experimen- 
tal situations in which the excitation typical period is 
large as compared to the magnetic diffusion time con- 
stant Tf ~ /i i 2 //Of, where pf stands for the flux flow 
resistivity and L is some typical length of the sample. 
Then, for current densities below J c , S vanishes as no 
electric field is generated in stationary conditions. On 
the other side, in the approximation of arbitrarily large 
flux flow resistivity, S would diverge if J > J c and, thus, 
J- n +i must be minimized constraining the current density 
to J < J c . We wish to remark that the definition of crit- 
ical current density may be done in a very general sense. 
In fact, one can postulate vertical E(J) relations for def- 
inite directions of space as it is the case of the DCSM, 
in which one uses parallel and perpendicular projections 
with respect to the local magnetic field. One could even 
dictate that huge dissipation occurs whenever the current 
density vector J lies outside some allowed region A, gen- 
crating an almost instantaneous change of the magnetic 
profile. 

In the light of the previous discussion, the evolution- 
ary critical state profiles can be obtained either by using 
Maxwell equations and a vertical E(J) law or the prin- 
ciple: 

In a type-II superconducting sample fl with an initial 
field profile H n (x), and under a small change of the exter- 
nal drive, the new profile H n+ i minimizes the functional 

C[H n+1 (x)} = ^ JjH n+1 - H n \ 2 , (4) 

with the boundary conditions imposed by the external 
source, and the constraint V x i? n +i € A(H n +\,x). 

In this work wc will concentrate on two particular cases 
which can be justified in terms of the underlying physi- 
cal mechanisms. First, we will analyze the isotropic case 



\J\ < Jc (i-e.: A is a disk). Then the variational formu- 
lation of the DCSM (Jn < J c u, J± < J c ±) will be given. 
In this case A is an oriented rectangle. 

The next section is devoted to some aspects on the 
mathematical formalism of the optimal control theory. 
We will show that this method provides a very convenient 
framework for applying the variational principle stated 
above. 



III. HAMILTONIAN EQUATIONS 

From the mathematical point of view, the problem of 
minimizing the action integral C[H n +i(x)] [Eq.(^)] with 
H n (x) given, and U n +\ fulfilling the differential equation 

V x H n+1 = J e A c M 3 (5) 

is a problem of variational calculus with nonholonomic 
constraints (constraints on the derivatives of the state 
variables). Being the constraint set A a compact region 
with boundary, we must take into account the possibility 
that the minimum is either reached in an interior point 
or in a boundary point. The machinery generalizing the 
classical variational calculus to sets with boundaries is 
the maximum -firinciple of optimal control, introduced 
by PontryaginEj OC is a standard theory in engineering, 
where dynamical systems are actuated by external con- 
trols (limited to some maximal values) in order to get a 
desired behaviour. The given system of differential equa- 
tions, including the effect of external actions is named 
the control system, and the integrand of the minimizing 
functional is the cost or performance function. The math- 
ematical problem here is similar, with Ampere's law as 
control system and performance function ^\H n+ i — H n \ 2 , 
although of course the bounded variable J is not an ex- 
ternal control. As we will see later on, the algebraic con- 
dition of maximality embraces the cases of interior and 
boundary optimal points, therefore including the classical 
equations of variational calculus for unconstrained prob- 
lems, and the modified equations for minimal solutiopa 
in the boundary. In addition to the original reference,E2l 
the interested reader is invited to review Refs. |23|,^4] for 
a more comprehensive and topical statement of the OC 
machinery. 

In this paper, in order to simplify the presentation, we 
fix the geometry of the sample to be an infinite slab, so 
that the control system (Ampere's law) becomes by sym- 
metry considerations a simple system of ordinary differ- 
ential equations. Other geometries, as the cylinder, give 
way to different control systems, and the case of a finite 
sample generates a more involved system of partial dif- 
ferential equations. For the slab geometry, with applied 
magnetic field parallel to the faces we have 

^p±l = f { H n+1 ,u,x) feA 
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as control system. Hereafter, we take the X axis per- 
pendicular to the slab faces and the origin of coordi- 
nates at the midplane. By construction, the vector 
/ = (0, f y , f z ) is orthogonal to the physical variable J, 
fy = Jz, fz = — Jy Notice that, despite the applied ro- 
tation, we use the same notation for the allowed control 
set A. Note also that the function / has dependence on 
the local magnetic field (this allows to include the usual 
models J C (\H\), and possible anisotropy), on the position 
(for potentially inhomogeneous materials) and on some 
independent coordinates u, the control variables in OC 
notation, parameterizing the region A. 

The first step in the theory of constrained variational 
calculus is the definition of a Hamiltonian density, con- 
taining the performance function, and momenta variables 
that play the role of Lagrange multipliers for the con- 
straints. 

H(H n+ i,u,p,x) = p- f- ^\H n +i - H n (x)\ 2 . (6) 

Recall that here H n (x) is a given profile. It is important 
to notice that H. is not the usual Hamiltonian function 
of Classical Mechanics in phase space; here it depends 
on the usual state H n+ i and momenta p variables, but 
additionally on the control u. Therefore, the associated 
equations are not only the Hamiltonian differential equa- 
tions but also an extra algebraic condition of maximality, 
in order to determine the extra variables u. Denoting by 
H* +1 (x), p*(x) and u*(x) the optimal solution functions 
(i.e. minimizing C and satisfying the control system), the 
OC equations are 



(7) 



-(H* +1 ,'0*,x) , (8) 



dH* +1 8U r , B , ^ , 

the adjoint equations for the momenta 

dp* _ dU 
~dx~~~dH n+1 ~ 

H* +l -H n {x)-p*-- 

and the algebraic condition of maximality 
H(H,u*,p,x)>H(H,u,p,x) V f(H,u,x) e A . (9) 

H and p are fixed in this last condition. In general, 
this allows to find a relation u*(H,p,x). When this re- 
lation is replaced in the former Hamiltonian equations 
[Eqs.(|7]) and (g)], a well posed system of ordinary differ- 
ential equations appears. 

For the class of Hamiltonian TC described above, the 
algebraic condition of maximality is fulfilled for a vector 
/ with maximum projection over the momentum p. As 
a simple exercise, the reader is invited to check that this 
rule produces the Bean CSM for one dimensional prob- 
lems (parallel vortices). Thus, on choosing H = H z z 



one has / = — J y . The restriction on the current den- 
sity reads J y — uJ c with |u| < 1. Then, Eq.(^) gives 
U* = sgn(p) and this leads to / = ± J c in agreement with 
Bean's prescription. Remarkably, for the case of multidi- 
mensional systems (non-parallel vortices) the maximal- 
ity condition will give the critical current vector both in 
modulus and direction. Note that J is determined dy- 
namically, through the evolution of the momentum vec- 
tor. The specific details on the distribution rule for the 
components of J c depend on the control space A and will 
be discussed in the next section. 

As regards the boundary conditions that must be used 
for obtaining the solution profiles H* +1 (x) from Eqs.(0) 
and (g), several considerations must be made. Firstly, in 
the absence of demagnetizing effects, H * +1 is determined 
on the faces of the slab by continuity of the external ap- 
plied field. For the remaining boundary conditions, two 
typical situations appear. In the first case, the modi- 
fied penetrating profile H* +1 (x) equals the former pro- 
file H n (x) before reaching the centre of the slab. This 
gives place to an (unknown in advance) point x* such 
that H* +1 (x) — H n (x) V x < x*. The free boundary 
condition Ti.(x*) — applies in this case, and allows to 
determine x* . In the second case, the new profile never 
meets the former one, and the value of H* +1 (0) is un- 
known; then, the corresponding transversality condition 
p(0) = completes then the number of required bound- 
ary conditions. 



IV. FIELD PENETRATION PROFILES 

The Hamiltonian formalism developed above can be 
expeditiously applied to calculate the field penetration 
profiles for magnetization processes in which non-parallel 
vortex configurations are enforced. As the applied field 
is assumed to be the same on both sides of the slab under 
consideration, we can restrict to the interval < x < a. 
By virtue of the symmetry, the same behavior appears 
for — a < x < 0. 



A. Isotropic model 

First, we will derive some results for the field pen- 
etration process in the so-called isotropic hypothesis: 
| J\ < J c . In order to show the capabilities of the model, 
a field dependence. J C (H) will be allowed. In particular, 
the Kim's modelcS expression J C (H) — J Co /(l + H/H ) 
will be used. Recall that the microstructure dependent 
parameters J C0 ,H are included. For convenience, the 
following dimcnsionless units are introduced: x is given 
in units of a, H in units of Hq, and J in units of Ho/a. 
Then, the statement of Ampere's law, together with the 
critical current restriction read 
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dHn+ 

dx 



f(H n+ i,u,x) 



ftu 



l + |-ff„ + l| 



l«l < 1 



Above we have introduced the dimensionless constant 
ft = J Ca a/Ho and the so-called control variable u, which 
is a vector within the unit disk D. This means that the 
current density belongs to a disk of variable radius, de- 
pending on the local magnetic field. Thus, following the 
OC terminology, we have the control equations for the 
state variables H n+ i(x). 

Next, we require the minimization of the functional 
C[H n +i(x)] constrained by the state equations. Follow- 
ing the OC machinery introduced in Sec.|l|, the maxi- 
mization of the associated Hamiltonian TL [Eq.(j^)] with 
respect to the control variable u leads to the condi- 
tion u* = p*/p*. This has the physical counterpart 
|J| = J C (H), i.e.: within the isotropic model, the maxi- 
mum allowed current density modulus J c is carried within 
those regions which have been affected by the perturba- 
tion. Then, the Hamiltonian equations that provide the 
field profile are 



dH, 



u+I.J = Pi_ 

p* 



ft 



dx 



dx 



n+l,i 



ftp*H* 



(10a) 
(10b) 



Eventually, appropriate boundary conditions must be 
supplied to solve this set of equations. As it was discussed 
in Sec. Ill these conditions are given by the external drive 
values, together with the field penetration regime. If the 
new profile matches the old one at a point < x* < 1 
(partial penetration) one uses H* +1 (x*) = H n (x*) and 
7i(x*) — 0. Within the full penetration regime, how- 
ever, one uses the momenta transversality conditions 
p*(0) = 0. 

In order to illustrate the bundle of phenomena which 
can be expected in crossed field measurements, we have 
calculated the penetration profiles for a zero field cooled 
sample to which a constant excitation H z s is applied, 
followed by cycling stages of the other field component 
at the surface H y s- We want to emphasize that, in or- 
der to approximate the continuum evolution, small incre- 
ments of H ys ha ve be en applied in the iterative solution 
of Eqs.(lOa) and ( |l0b| ). However, for clarity, only a selec- 
tion of representative field and current density curves are 
depicted. The value of ft is set to unity for definiteness 
hereafter. 



1. Initial magnetization process 

Fig.|l| displays the penetration profiles for the initial 
magnetization process that is induced by increasing the 
surface field H y s- We notice that the conventional H z (x) 
profile (which was obtained by the standard CSM) is 



pushed towards the center as H y s increases. In phys- 
ical terms, a current flow J z is required for shielding 
the new field component H y . Owing to the constraint 
| J| = J C (H), this results in a reduction of J y (the slope 
of H z ) and H z develops a subcritical behavior. Just for 
clarity, we have only plotted advancing fronts until the 
centre of the sample is reached. Ongoing profiles will be 
shown in the next subsections. 



2. Low field hysteresis 

Here, we analyze the hysteresis effects which can be 
observed when the surface field is cycled. First, the max- 
imum applied field (—H m < H y s < H m ) is chosen so as 
to keep a zero field value at the midplane. We define the 
penetration field as the value of H y s for which the 
front reaches the point x = 0. This is somehow equivalent 
to the so-called partial penetration regime in the stan- 
dard CSM. However, notice that now depends on the 
previously set value of H z s- FigJ| displays the results. 
It can be observed that H y (x) basically shows a con- 
ventional critical state behavior, whereas H z (x) develops 
a prominence. This shape relates to the corresponding 
slope reduction near the surface (x = 1) and the subse- 
quent reentry in the unperturbed profile. Physically, the 
distribution rule for the current density vector enhances 
the J z component near the surface in order to shield the 
change in H y (i.e.: reducing J \H n+ i — H n \ 2 dx). Then, 
J y diminishes and, therefore H z (x) flattens. An inter- 
nal slope increase is obviously required for the previous 
profile H n _ z (x) to be reached at some point < x* < 1. 
Notice the peak structure that appears in J y (x) and the 
corresponding sign change in the J z component. The 
point of vanishing J z determines the maximum of the 
peak structure in J yi i.e: the full current flow is parallel 
to the Y axis. 



3. Intermediate field hysteresis 

When the maximum applied field H m is allowed to 
increase beyond new phenomena appear. Fig.|| dis- 
plays some outstanding features, which are sketched here: 
(i) it is apparent that the H y (x) curves do not show a 
conventional hysteresis cycle structure. Notice that the 
profiles no.l and no. 11 that correspond to H y g = H m 
do not fit. (ii) Contrary to the case of low field cy- 
cles, the profiles H z (x) do not repeat values for the as- 
cending/descending branches. When the sample centre 
is reached (profile no. 4) irreversible flux entry happens 
(no. 5 and no. 6). Further changes (curves no. 7 to no. 11) 
occur over a reduced shielding curve. Finally, one can see 
that the aforementioned features are clearly translated to 
the current density vector J(x). 
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4- High field hysteresis 

In this part we show that, within the isotropic model, 
a transverse field H y s may effectively collapse the lon- 
gitudinal diamagnetic profile H z (x). Fig.|| displays a 
magnetization process in which H y g is increased up to 
the value H m = 3 in our dimensionless units. Then, 
H y s is decreased down to H rn — —3. Further increase to 
H m = 3 is not shown just to avoid intricacy in the graph. 
However, the complete hysteresis cycle will be presented 
in the next section. Here one can see that H y (x) fol- 
lows a typical CSM high field evolution, whereas H z (x) 
irreversibly evolves to a flat profile, compatible with the 
boundary condition H z (l) = H z s- Simultaneously, the 
current density related to H z becomes zero. Thus, a con- 
stant trapped field H z g is generated within the sample. 
Any subsequent change in H y s will not affect the trapped 
field and the associated current density component. 

The described behavior can be well understood in 
terms of the minimization principle. The process tends 
to minimize 



C 



+l,y 



(■ffn+l,z — Hn,z) 



It is apparent that if the condition H n+ i, z (x) = H n ^ z (x) is 
fulfilled, C cannot be smaller than the minimum value of 



1/2 f(H n +i t y — H nt y) 2 , which appears for the full current 
flow dedicated to H v . Thus, one has J z = ±J C and J v = 



l+l. y "n,y 

L V . Thus, one has J. 
0, and both conditions hold thereafter. 



B. Double critical state model 

Below we will show that our approach reproduces the 
results derived by the standard DCSM when the appro- 
priate control space for the current density is chosen. For 
this purpose we focus on Ref. [l2]in which the DCSM was 
used for investigating a wealth of phenomena in rotating 
field experiments. In order to ease comparison, we shall 
adopt the notation therein. An infinite slab with surfaces 
at x = and x = 2X m was considered with the magnetic 
field contained in the YZ plane. 

The starting point will be the restatement of the func- 
tional C in appropriate generalized coordinates 

C[H n+1 (x)} = - I - 2H D H n+1 cos (a n+ i - a n ) . 

1 Jn 

Here H, a stand for the field modulus and the angle with 
respect to the Y axis. H n ,a n correspond to the given 
previous profile and one must solve for H n+ i, a n +i. Now, 
C must be minimized with the constraints 



I dH n+ \ 
1 dx 1 



< J, 



c_L 



I da n+1 
1 dx 



< k c \\ , 



where J c ± and fc c || are the constants characterizing the 
flux pinning and cutting thresholds. Following Rcf. n3 



fc c || relates to the parallel critical current density by fc c || = 
J c \\/H. 

On defining the characteristic length scale A = l/fc c ||, 
and using H in units of J c j_A and x in units of A, one 
obtains the dimensionless expressions 



dH n+ i 


Jx 


dx 


JcX 


da n+ i 


= k L 


dx 


k c \\ 



= u h 



\u h \ < 1 



In other words, the control set A is a square in this units. 

Next, we introduce the associated Hamiltonian [Eq. (||) 
in Sec.|l| 

H = phUh +p a u a - - - 2H D H n+1 cos (a n+ i - a n )] . 

Then, if one applies the Pontryagin's maximum principle 

max7i = max(phUh + p a u a ) 
ueA SeA 

is determined again as a vector leaning on the 
boundary of A. This has been illustrated in Fig.|| and 
relies on the fact that p ■ u is maximum for the largest 
projection of u on p. If (j>h,Pa) are non-vanishing, one 
gets {u* hl u* a ) = [sga{ph), sgn(p Q )], which can be iden- 
tified as a CT-zone within the DCSM framework. On 
the other hand, ph ^ 0,p a = in an open interval will 
give place to a T-zone in which u* h = sgn(p^) and u* 
must be determined on the basis of other arguments. Fi- 
nally, ph — 0,p a ^ corresponds to a C-zone, in which 
it* = sgn(p a ) and Uh* is yet undetermined. 

Thus, in the case under consideration (i.e.: field inde- 
pendent J c ±, fc c |i) the Hamiltonian equations are 



dH n+ i 
dx 

da n +i 
dx 

dph 
dx 

dp a 

dx 



= sgn(ph) or undetermined (H a ) 

= sgn(p Q ) or undetermined (lib) 

= H n+1 - H n cos (otn+i - oc a ) (11c) 

= H n+1 H n sm(a n+ i - a„). (lid) 



These equations can be straightforwardly solved to pro- 
duce the field penetration profiles. Different expressions 
arise related to the mentioned zone structure. 



1. CT zone 

If one has ph ^ 0,p a =/= 0, the field penetration is given 
by the differential equations 



dH n+1 
dx 

da n +i 
dx 



= sgn(p/,) =>■ H n+1 = sgn(p h )x + H n+1 (0) 
= sgn(p a ) ^ a n+ i = sgn(p h )x + a n+ i(0) . 



G 



Above, we have incorporated surface boundary condi- 
tions, just for definiteness. 

2. T zone 

In this case one has 

H n+1 = sgn(p h )x + H n+1 (0) 

and a n +\(x) is determined by the condition dp a /dx = 0, 
which guarantees that p a keeps the zero value in this 
region. Then 

H n+1 (x)H a (x) sin [a n+1 (x) - a n (x)] = 

3. C zone 

In this case one has 

a n+1 = sgn(p h )x + a n+ i(0) 

and H n+ i(x) is determined by the condition dph/dx = 0, 
which gives 

H n+1 (x) = H n (x) cos [a n+ i(x) - a n (x)} 



and defines the beginning of a C-zone. The equation de- 
termining x Vtn+ i has been obtained within the approxi- 
mation 8a — > cos™ (8a) — > 1. This condition yields 
the continuum limit of our discretized approach. No- 
tice that our expression matches the one obtained by the 
conventional DCSM. Finally, the C-zone stretches from 
to 

x c ,n+ii the point where a n +i vanishes. Then, 
the held modulus profile becomes 

H n+1 (x) = H s cos (a s ,„+i - x) . (12) 

The associated momenta PhiVa may be straightforwardly 
obtained by integration of their first order differential 
equations. In Fig.^j we have plotted the first step pro- 
files Hx(x), a±(x), ph(x),p a (x). Several points have been 
marked on the curves Ph{x) and p a (x), that correspond 
to the sequence of values of the control variable u stand- 
ing out in Figj|. Notice that (tih, u a ) always takes values 
on the boundary of the control space. Within a CT- 
zone u remains on a vertex of the square. As soon as 
a component of p vanishes (point 3 in our example) the 
associated component of u must be determined by the 
condition of zero derivative. In our case, this produces a 
jump to the vector u^, which defines the beginning of the 
C-zone (U3, M4, U5, . . .). This jump in the control vari- 
able is known as bang-bang phenomenon among mathe- 
maticians. Eventually, if the center of the sample is not 
reached by the variations, a new jump in the current 
density gives way to the O-zone (uh — 0, u a = 0). 



4. Example: nonmagnetic initial state 



V. MAGNETIZATION LOOPS 



Several initial magnetic configurations were examined 
in Ref. [O] within the conventional DCSM formulation. 
For illustration, the nonmagnetic initial state, which cor- 
responds to the so-called field cooled experiment is ana- 
lyzed here in the framework of our theory. 

We assume initial uniform flux density within the 
superconductor Hq[x) = H s ,a$(x) = 0. Then, 
consider the quasisteady evolution towards the state 
H n +i(x), a n+ i(x) with boundary conditions at the sur- 
face H n+ i(0) = H s and a n +i(0) = (n + l)8a = a s ,„+i. 
The minimization of C for the n + 1 time layer requires 
an initial CT-zone given by 

H n+1 (x) = H s - x ; p h < 
a n +i(x) = a StTl+ i - x ; p a < . 

The CT-zone spreads up to the point x v ^ n +\ where the 
conditions p^ — 0, dp^/dx = are reached. This point 
may be determined by iteration of the governing Hamilto- 
nian equations (Eqs. (fil|)) over the previous steps. Such 
process leads to 

H s - x v<n+ i - H n (x tun+ i) cos (8a) =0 => 
%v,n+i = H s [1 - cos (a s>n+ i - x v , n+ i)] 



Though some experiments have been devised for the 
direct measurement of the field or current penetration 
profiles, the common practice is to record average values. 
Thus, the vast majority of data on vortex matter in type 
II superconductors consist of magnetic moment measure- 
ments. Recall that the macroscopic magnetization of the 
sample is defined as M = (H(x)) — Hs, where Hs stands 
for the uniform applied field. 

In this section we concentrate on the magnetization 
loops which can be derived from the previously treated 
field penetration profiles. It will be shown that important 
experimental results can be reproduced. New observ- 
able features are also predicted, which may be used as 
a test for the physically meaningful critical current con- 
trol space. All the results presented below correspond to 
the isotropic hypothesis, which has been mainly devel- 
oped in this paper. The comparison to the already well 
established DCSM and a discussion on the selection of 
the control space for the current density are considered 
in the next section. 

First, we analyze the influence of the maximum ap- 
plied field H m on the magnetization loop. Fig.0 col- 
lects the results corresponding to the profiles analyzed 
in SubSec.IV A: low field (H m < H^), intermediate field 
(H m > H* n ), and high field (H m > H* n ). We have 
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depicted both M y and M z for the zero field cooled pro- 
cess: (i) Apply H z s, (ii) increase H y s from zero to H m , 
(iii) decrease H y s to —H m , and (iv) increase H y s to H m 
again. 

My(H y s) roughly displays a standard behavior, but 
shows a noticeable peculiarity for the intermediate field 
loop. The loop fails to close after a complete field cy- 
cle. This feature was not noticed for the low/high field 
loops, at least to the numerical precision of our calcu- 
lations. The M z (H y s) curves show a quite less conven- 
tional behavior as compared to one dimensional measure- 
ments in which a single field component appears. In ad- 
dition, qualitatively different tendencies can be observed 
in terms of the value H m . Thus, low field cycles produce 
a butterfly loop; intermediate fields produce a descending 
behavior for the transverse magnetic moment, and high 
fields produce the irreversible collapse of magnetization. 

As a final observation, we describe an outstanding fea- 
ture that was observed in our simulations. A second cycle 
in the applied field H y s succeeds to close the hysteresis 
loop M y (H y s) in the intermediate field region. Mean- 
while, M z keeps the descending trend, though seeming 
to stabilize around a non-vanishing value. This has been 
depicted in Fig.[| 

We want to emphasize that our intermediate field sim- 
ulations reproduce the experimental observations in Ref. 
^ and in some related papers cited therein. Thus, the 
isotropic hypothesis sets a simple model for those crossed 
field measurements. Furthermore, we suggest to extend 
such measurements as a test of the theory, i.e.: perform 
low/high field loops and a second cycle in order to check 
the predicted features. 

VI. DISCUSSION AND CONCLUSIONS 

We have introduced a variational formulation that gen- 
eralizes Bean's critical state model for hard superconduc- 
tors, while keeping its conceptual simplicity. In the spirit 
of the original theory, we have also avoided the explicit 
use of the electric field E. However, the time evolution of 
the magnetic field profiles H (x, t) can be obtained by re- 
placing E with an appropriate restriction on the current 
density J flowing within the sample. Mathematically, a 
wide range of possibilities are allowed as we pose it in the 
form J € A(H,x) with A a bounded set. Upon the se- 
lection of a definite set A, our theory may be considered 
as the background for a variety of critical state models. 
In particular, we have analyzed the cases of choosing ei- 
ther a circle or a rectangle (isotropic model and DCSM 
respectively). 

The variational quantity which is minimized under the 
aforementioned restriction has got a clear physical inter- 
pretation and quite general validity. We propose that 
the time evolution is governed by the balance between 
a screening term (inertia) and an entropy production 
term. This relies on basic principles of non-equilibrium 



thermodynamics.^ Within the usual assumption of in- 
stant response (vertical E(J) graph) the singular behav- 
ior of S = E ■ J/T is treated by replacing the associated 
term with the condition J e A in the minimization prin- 
ciple. This kind of statement fits the mathematical the- 
ory of optimal control,EJ that has been applied in this pa- 
per. The method is appropriate for the numerical treat- 
ment of realistic problems in which the magnetization is 
a vector quantity. Such situation occurs in experimental 
setups designed for studying the interaction of twisted 
vortices and in the field of technological applications of 
hard type-II superconductors. 

We want to emphasize that the selection of the permit- 
ted set A for the current density should be either justified 
in terms of fundamental theories or bounded by experi- 
mental observations. The scope of our paper is to pro- 
vide a suitable mathematical treatment for general CSM 
problems and a physical background for the method. A 
well posed inversion theory, allowing to obtain the set A 
from experiments is still an open question. 

However, some conclusions can be drawn from the 
comparison of simulations with different sets A. On 
the one side, some experimental observations are quite 
independent of the region A in use. This is the case 
of the field penetration profiles in rotating field exper- 
iments. A very similar behavior is_predicted either by 
the isotropic or DCSM hypothesis.E£l On the other side, 
one can find outstanding differences between these two 
models. Firstly, the former predicts the criticality con- 
dition J = J c for those regions which have been affected 
by an electric field. Within the DCSM framework, how- 
ever, one can obtain critical regions (CT zones) where 
J = J c = \/(Jc\\ + Jq±) an d subcritical regions (C or T 
zones) where J < J c . Microscopic probes can be used to 
investigate this property. For instance, Ref. |26| reports 
neutron polarization analysis of the magnetic field pro- 
files within ceramic YBCO samples subjected to rotation 
stages. The isotropic hypothesis seems to be favoured for 
this Josephson medium. On the macroscopic side, one 
can design adequate experiments which check the com- 
patibility with specific selections of A. The investigation 
of the so-called magnetization collapseu is a promising 
possibility. In fact, recent experimental rcsiiits on melt- 
textured and single crystal YBCO samplesuEj are nicely 
described within the isotropic model. The authors of Ref. 

have studied the transverse magnetic moment M z sup- 
pression by cycling the longitudinal magnetic field H y s- 
As the most important observed feature they emphasize 
the symmetric decrease of M z for the diamagnetic and 
paramagnetic initial states of the sample. This is shown 
to be in clear contradiction with the DCSM hypothesis 
and qualitatively interpreted by their two-velocity hydro- 
dynamic model. In Fig.^jwe show that our general critical 
state theory gives a quite satisfactory explanation of the 
experimental facts. We get a symmetric suppression of 
M z . The field penetration profiles show that the slope 
decrease in H z (x) near the surface, due to the enhance- 
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ment of the current density component which shields the 
change in H y , occurs in a quite similar fashion, regard- 
less the initial magnetic state. Additionally, we want 
to notice that the behavior of M z upon a second cycle 
in H y s which is observed in our simulations (see Fig j^) 
gives further evidence of the validity of the model. The 
coincidence with the experimental data depicted in Fig. 5 
of Ref. [2£] is remarkable. 

Finally, we will comment on some future extensions 
of our work. Here we have used a time-discretization 
scheme, and a slab geometry, in order to produce an or- 
dinary differential equation statement of our principle. 
This simplifies the theory, but should not be considered 
as a limitation. Optimal control for systems governed 
by partial differential equations may be applied for more 
general cases. Thus, one could afford a continuum ap- 
proach to the problem and also include the sample finite 
size effects if necessary. Other issues as the use of differ- 
ent physically meaningful control sets A (f.L: elliptic, so 
as to incorporate material anisotropy), spatial inhomo- 
geneity or surface barrier effects are also suggested. 
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FIG. 1. Field and current density penetration profiles for 
a zero field cooled slab. The dots indicate the initial process 
in which H z was raised to the value H z s = 0.6 at the surface. 
The successive steps increasing H y s are labelled from 1 to 
5. Continuous lines are used for H z (x) and the associated 
current density —J y (x), while dashed curves represent H y (x) 
and Jz{x). Dimensionless units are used as defined in the 
text. 



FIG. 2. Field and current density penetration profiles 
for a small amplitude cycle in H y s, successive to the initial 
magnetization process described in Figjl| The descending 
branch curves are labelled 1 to 7 (continuous lines), whereas 
ascending curves are labelled 8 to 14 (dash-dotted curves). 
The Hz(x) and —J y (x) profiles nearly repeat values for the 
descending/ascending processes as one can see in the graph. 
Some labels have been skipped to avoid confusion. 







FIG. 3. Field and current density penetration profiles for 
an intermediate amplitude cycle in H y s, successive to the 
initial magnetization process described in Fig.[j]. The curves 
labelled 1 to 6 (continuous lines) stand for the descending 
process, while the curves 7 to 11 (dash-dotted) stand for the 
ascending part. The label no. 7 in the family of curves H z (x) 
has been avoided for clarity. 



FIG. 4. Field and current density penetration profiles 
corresponding to a high amplitude cycle in H y s, successive to 
the initial magnetization process described in Figjl|. Here, we 
have labelled 1 to 7 the profiles for the initial increase to the 
high amplitude {H y s = 3). The descending branch curves are 
labelled 8 to 14. In order to avoid confusion, the family of 
ascending branch curves of the cycle are not displayed. 



FIG. 8. First and second magnetization cycles for the 
intermediate field loop considered in FigQ 



FIG. 9. Dependence of the transverse magnetization M z 
on the first cycle of H y s (0 — * H m — > — H m — > H m ) for 
the para- and diamagnetic initial states in H z (x). The upper 
panel shows the suppresion of M z in the paramagnetic case 
(sequence O'A'B'C) and the lower panel reflects the same 
phenomenon for the diamagnetic case (sequence OABC). The 
insets illustrate some selected penetration profiles of H z (x), 
corresponding to the points labelled on the M z curves. 



FIG. 5. Diagram of the square control space A on purpose 
for the DCSM oriented according to the local magnetic field. 
For illustration, we plot the optimal control solution u which 
maximizes p ■ u in a generic situation. Some specific values 
pi , P2 , P3 , ■ ■ ■ and the corresponding solutions ui, u<2, Us, ■ . . are 
also marked. The sequence is extracted from the example 
depicted in FigH 



FIG. 6. Normalized field modulus and angle penetration 
profiles (H(x)/Hs,a(x)/cts) obtained by the application of 
our variational principle to the DCSM control space. A field 
cooled slab was assumed, subjected to further rotation of 
the surface field. Also depicted are the associated momenta 
(ph{x),p a (x)). The points 1 to 5 labelled on the momenta 
curves correspond to the sequence of control values shown in 
Fig.^. Following the notation introduced in Ref. we have 
marked the CT, C and O zones, as well as the penetration 
points x v and x c . The distance x is given in units of the slab 
half- width x m . Contrary to our choice in the rest of figures, 
in this case x — corresponds to the surface of the slab. 



FIG. 7. Evolution of the magnetization components in a 
simulated crossed field experiment for a zero field cooled su- 
perconducting slab. Subsequent to the application of a con- 
stant surface field H z s, the other component H y s was cycled 
(0 — * H m — > — H m — ► Hm) for three different values of H m . 
Continuous lines stand for low fields, continuous-dotted for 
intermediate fields and dashed lines represent the high field 
cycle. 
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